IMPLEMENTING THE ASYMPTOTICALLY FAST VERSION OF THE 
ELLIPTIC CURVE PRIMALITY PROVING ALGORITHM 


F. MORAIN 


ABSTRACT. The elliptic curve primality proving (ECPP) algorithm is one of the current fastest 
practical algorithms for proving the primality of large numbers. Its running time currently cannot 
be proven rigorously, but heuristic arguments show that it should run in time O((log N)°) to prove 
the primality of N. An asymptotically fast version of it, attributed to J. O. Shallit, is expected to 
run in time O((log N)*). We describe this version in more details, leading to actual implementations 
able to handle numbers with several thousands of decimal digits. 


1. INTRODUCTION 


From the work of Agrawal, Kayal and Saxena [2], determining the primality of an integer N can 
be done in proven deterministic polynomial time O((log N)!°5). More recently, H.-W. Lenstra, Jr. 
and C. Pomerance have announced a version in O((log N)®) (see [31]). Building on the work of 
P. Berrizbeitia [7], and reusing classical cyclotomic ideas that originated in the Jacobi sums test 
(1, 11], D. Bernstein [6] and P. Mihailescu & R. Mocenigo [34], independently, have given improved 
probabilistic versions with a claim of proven complexity of O((log N)*). For more on primality 
before AKS, we refer the reader to [14] (see also [36]). For the recent developments, see [5]. 

All the known versions of the AKS algorithm are for the time being too slow to prove the primality 
of large explicit numbers. On the other hand, the elliptic curve primality proving algorithm [3] has 
been used for years to prove the primality of ever larger numbers*. It is sketched in [30] that 


Conjecture 1.1. The cost of ECPP is O((log N)°). 


The same article describes an asymptotically fast version of ECPP, attributed to J. O. Shallit 
and that we will call fastECPP: 


Conjecture 1.2. The cost of fastECPP is O((log N)*). 


The aim of the present paper is to describe fastECPP, give a heuristic analysis of it and describe 
its implementation. 

Section 2 collects some well-known facts on imaginary quadratic fields, that can be found for 
instance in [13]. Section 3 presents the basic ECPP algorithm and analyzes it. In Section 4, the 
fast version is described and its complexity estimated. Section 5 explains the implementation and 
Section 6 gives some actual timings on large numbers. All estimates make use of the O notation. 
Remember that O(f) means O(f (log f)O°™). 
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2. QUADRATIC FIELDS 


A discriminant —D < 0 is said to be fundamental if and only if D is free of odd square prime 
factors, and moreover D = 3 mod 4 or D mod 16 € {4,8}. The quantity 


D(X) = #{D < X,—D is fundamental} 


is easily seen to be O(X). 
A fundamental discriminant may be written uniquely (modulo permutations) as: 


t 
-D=||¢ 
4=1 


where q* is either —4 or +8, or g¥ = (—1/q;)q; for any odd prime q;. (We denote by (4) the Legendre 
a a DP. 


symbol.) The number of genera is g(—D) = 2°"! and Gauss proved that this number divides the 
class number h(—D) of K = Q(./—D). Moreover, Landau (see e.g., [37, Thm 4.4, p. 143]) proved 
that h = O(D'/? log D). 

The rational prime p is the norm of an integer in K, or equivalently, 4p = U? + DV? in rational 
integers U and V if and only if the ideal (p) splits completely in the Hilbert class field of K, 
denoted K 7, an extension of degree h(—D) of K. The probability that a prime p splits in K 7 is 
1/(2h(—D)). ; 

Using Gauss’s theory of genera of forms, it is known that if (=) = 1 for all 2 (equivalently, (p) 
splits in the genus field of K), then the probability of (p) splitting in K y is g(—D)/h(—D). 


3. THE BASIC ECPP ALGORITHM 


We present a rough sketch of the ECPP algorithm, enough for us to estimate its complexity. We 
do not insist on what happens if one of the steps fails, revealing the compositeness of N. More 
details can be found in [3]. Without loss of generality, N is supposed to be greater than 3. 


3.1. Elliptic curves over Z/NZ. For us, an elliptic curve E modulo N is an equation Y? = 
X°+aX +b, a,b€ Z/NZ with gcd(4a? + 27b?, N) = 1 and we will let 


E(Z/NZ) = {(a:y:z) € P?(Z/NZ), y?z = 2? + axz? + bz*} U {Og = (0:1:0)} 


which ressembles the definition of an actual elliptic curve if N is prime, P?(Z/NZ) being the 
projective plane over Z/NZ. The important point here is that if p is a divisor of N, we can consider 
E as acurve E, over Z/pZ and reduce a point P € E(Z/NZ) to Py € Ep(Z/pZ). Moreover, we can 
define an operation on E(Z/NZ), called pseudo-addition, that adds two “points” P and Q with 
the usual chord-and-tangent law. This operation either yields a point R or a divisor of N if any is 
encountered when dividing. If R exists, then it has the property that Rp is the sum of P, and Q, 
on E, for all prime factors p of N. Note also that Og reduces to the ordinary point at infinity on 
Eg: 

We will need to exponentiate points in EF. This is best defined using the division polynomials 
(see for instance [4] for a lot of properties on these). Remember that over a field K there exists 
polynomials $,,(X,Y), Um(X,Y), wm(X,Y) such that, if P= (xX :Y:1), 

(1) [m]P = P+---+P = (bm(X,Y)tm(X,Y) : wm(X,Y) + ¥3,(X,Y)). 

m times 
All these polynomials can be computed via recurrence formulas and there is a O(log m) algorithm 
for this task (a variant of the usual binary method for exponentiating). A fundamental property is 


the following: 
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Theorem 3.1. Let K stands for an algebraic closure of K. Let m be an integer. If P = (x: y: 


1) € E(K) is such that [2)}P 4 Og, then P is of m torsion if and only if Um(x,y) = 0. 
For the sake of presenting the algorithm in a simplified setting, we prove (compare [27]): 


Proposition 3.1. Let N’ a prime satisfying (WN — 1)? <2N’ <(/N+1)?. Suppose that E is a 
curve modulo N, that P = (x: y: 1) € E(Z/NZ) is such that gcd(y, N) = 1, Wan (a, y) = 0 mod N 
but gcd(wyr(x,y),N) =1. Then N is prime. 


Proof: suppose that N is composite and that p < VN is one of its prime factors. Let us look at 
what happens modulo p. Since p{ y, Py is not a 2-torsion point on Ey. Since Wy (x, y) is invertible 
modulo p, then [N’](P,) 4 Oz, and therefore P, has order 2N’ on Ep. This is impossible, since 


IN! > VN = 1)? = p—1)? s (/p + 1)? > #E, by Hasse’s theorem. 


3.2. Presentation of the algorithm. We want to prove that N is prime. The algorithm runs as 
follows: 


[Step 1.] Repeat the following: Find an imaginary quadratic field K = Q(./—D) of discriminant 
—D, D> 0, such that 


(2) 4N =U? + DV? 


in rational integers U and V. For all solutions U of (2), compute m = N 4+ 1 — U; if one of these 
numbers is twice a probable prime N’, go to Step 2. 


Step 2.] Build an elliptic curve E over Q having complex multiplication by the ring of integers O x 
of K. 


Step 3.] Assuming N is prime, it splits completely in Ky. Reduce E modulo a prime divisor of 
(NV) in Ky, to get a curve modulo N. 

Step 4.] Find P = (a: y: 1), gcd(y, N) = 1 0n E such that Won (a, y) = 0, but ged(wyr (a, y), N) = 
1. If this cannot be done, then N is composite, otherwise, it is prime by Proposition 3.1. 


Step 5.] Set N = N’ and go back to Step 1. 


3.3. Analyzing ECPP. We will now analyze all steps of the above algorithm and give complexity 
estimates using the parameter L = log N. One basic unit of time will be the time needed to 
multiply two integers of size bounded by L, namely O(L!*“), where 0 < p< 1 (4 = 1 for ordinary 
multiplication, or any « > 0 for fast multiplication). 

Clearly, we go through Step 5 O(L) times for proving the primality of N. We consider all steps, 
one at a time, easier steps first. 


3.4. Analysis of Step 4. Finding a point P can be done by a simple algorithm that looks for the 
smallest x such that x? + ax + b is a non-zero square modulo N and then extracting a squareroot 
modulo N, for a cost of O(L?*“), using for instance the Tonelli-Shanks algorithm (or Lehmer’s 
algorithm when N — 1 is divisible by a large power of 2). Using the trick described in [3, §8.6.3], 
we do not need the modular square root, though... 

Computing yy (a, y) costs O(L?*“), and we need O(1) points on average, so this steps amounts 
for O(17"*), 
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3.5. Analyzing Step 2. The original version is to realize Ky/K via the computation of the 
minimal polynomial Hp(X) of the special values of the classical j-invariant at quadratic integers. 
More precisely, we can view the class group Cl(—D) of K as a set of h primitive reduced quadratic 
forms of discriminant —D. If (A, B,C) is such a form, with B? — 44C = —D, then 


Hp(X)= JJ (X= i((-B + V=D)/(24))). 
(A,B,C)€Cl(—D) 
It is known that Hp(X) € Z[X]. In [16, 15], it is argued that the logarithmic height of this 
polynomial (logarithm of the sup norm) is well approximated by the quantity DS(D) where: 


S(D)= Ss" 2. 
[A,B,C]eCl(—D) 
It can be shown that S(D) = O((log h)?). 
Evaluating the complex roots of Hp(X) and building this polynomial can be done in O(D) 
operations (see {15]). Note that this step does not require computations modulo N. 
Alternatively, we could use the method of [12, 8] for computing the class polynomial and get a 
proven running time of O(h?), but assuming a suitable Riemann hypothesis. 


3.6. Analyzing Step 3. To get E modulo N, we need a root of Hp(X) modulo N. This can 
be done with the Cantor-Zassenhaus algorithm (see [23, §14.3] for instance). Briefly, we split 
recursively Hp(X) by computing ged((X + a)(N~)/? — 1, Hp(X)) mod N for random a’s. 

Computing (X+a)(N—)/? mod (N, Hp(X)) costs O(LM(N, h)) where M(N, d) is the time needed 
to multiply two degree d polynomials modulo N. A gcd of two degree d polynomials costs 
M(N, d) log d (see [23, Ch. 11]). 

We expect a splitting after O(1) trials; after each splitting, we carry on with the factor of lowest 
degree. Since the degree is divided by at least 2 after each splitting, we get a root after O(log h) 
steps, but the overall cost is dominated by the first step, for a cost of 


O(M(N, h)(L + log h)). 
We can assume that M(N,d) = O(d'*”L'*¥) where again 0 < v < 1. 


3.7. Analysis of Step 1. This is the crucial step. Given D, testing whether (2) is satisfied involves 
the reduction of the ideal (N, roy=D) that lies above (NV) in K, where r? = —D mod (AN) (if N is 
prime... ). This requires the computation of V—D mod N as in Step 4, i.e., a O(L?*") time. Then 
it proceeds with Cornacchia’s algorithm, akin to a half gcd (see 5.2), for a cost of O(L'*¥). 

In the event that equation (2) is solvable, then we need to check that m = 2N’ and test N’ for 
probable primality, which costs again some O(L?*¥). 

The heuristic probability of m being of the given form is O(1/L). Though quite realistic, this 
estimate is out of reach of current techniques in analytical number theory. Using this heuristic, we 
expect to need O(L) discriminants D. Let us take all discriminants less than Dyax. Since N splits 
completely in Ky with probability 1/(2h(—D)), we expect to get a useful m provided 


1 
Se eas, 
eee 2h(—D) 
The left hand side dominates 
1 V Dyax 
Ved WeDo 
ae D log D 0g Umax 


so Dmax = O((Llog L)”) = O(L) should suffice. 
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Turning to complexity, Step 1 then solves O(L?) equations of type (2), followed by O(L) probable 
primality tests: 


O(?( Ete SLT y+ OF. Tere 3%, 
V—D mod N _ reduction probable primality 


which is dominated by the first cost, namely O(L4*“). 


3.8. Adding everything together. Taking D = O(L?) readily implies h = O(L). The cost. of 
Step 2 is O(L”), and that of Step 3 is O(L°+#*”), which dominates Step 4. All in all, we get that 
ECPP has heuristic complexity O(L**") for one step, and therefore O(L°+“) in totality. 


3.9. Remark. In practice, the dominant term of the complexity of Step 1 is O(n pL?*“) where np 
is the number of D’s for which we solve equation (2). Depending on implementation parameters 
and real size of N, this number np can be quite small. This gives a very small apparent complexity 
to ECPP, somewhere in between L? and L* and explains why ECPP seems so fast in practice (see 
for instance [22]). 


4. THE FAST VERSION OF ECPP 


4.1. Presentation. When dealing with large numbers, all the time is spent in finding the discrim- 
inants D, which means that a lot of squareroots modulo N must be computed. A first way to 
reduce the computations, alluded to in [3, §8.4.3], is to accumulate squareroots, and reuse them, 
at, the cost of some multiplications. For instance, if one has —3 and V5 = /—20//—4, then we 
can build /—15, etc. 

A better way that leads to the fast version consists in computing a basis of small squareroots 
and build discriminants from this basis. From the analysis carried out above, we need O(L?) 
discriminants to find a good one. The basic version finds them by using all discriminants that are 
of size O(L?). As opposed to this, one can build those discriminants as —D = (—p)(q), where p 
and q are taken from a pool of size O(L) primes. 

More formally, we replace Step 1. by Step 1’. as follows: 


[Step 1’.] 
1.1. Find the r = O(L) smallest primes q such that (5) = 1, yielding O = {q7,93,..., ar}. 
1.2. Compute all /q¢* mod N for q* € Q. 
1.3. For all pairs (qj, ,q;,) of Q for which qj qj, = —D < 0, solve equation (2). 


7 The cost of this new Step 1 is that of computing r = O(L) squareroots modulo N, for a cost of 
O(L- L?*#). Then, we still have O(L?) reductions. The new overall cost of this phase decreases 
now to: 


+ : 2+ (72, l+p 4+ . 2+ — O(7 3+ 
OL. L+H )4+ O(n? DH )+O(L. OL ) = O(L3+#), 
squareroots reduction probable primality 


Note here how the complexity decomposes as 3 = 1+ 2 or 2+ 1 depending on the sub-algorithms. 
Putting everything together, we see that this time the dominant step in terms of complexity is 
that of Step 3 and therefore we end up with a total cost of O(L**#*”) for this variant of ECPP 


and O(L*) asymptotically. 


4.2. Remarks. 


4.2.1. Complexity issues. We can reduce further the number of square root computations by using 
all subsets of Q and not only pairs of elements. This would call for r = O(log log N), since then 
2” = L? could be reached. Though useful in practice, this phase no longer dominates the cost of 
the algorithm. 

Moreover, we can see that several phases of fastECPP have cost O(L?), which means that we 


would have to fight hard to decrease the overall complexity below O(L*). 


4.2.2. A note on discriminants. Note that we use fundamental discriminants only, as non funda- 
mental discriminants lead to curves that do not bring anything new compared to fundamental ones. 
Indeed, if D = f?D, with D fundamental, then there is a curve having CM by the order of discrim- 
inant D. Writing 4N = U2? + Df?V?, its cardinality is N + 1—U, the same as the corresponding 
curve associated to D. 


4.2.3. A note on class numbers. As soon as we use composite discriminants —D of the form qj. qj, 
Gauss’s theorem tells us that the class number h(—D) is even. This could introduce a bias in our 
estimation, but we conjecture that the effect is not important. 


5. IMPLEMENTATION 


5.1. Computing class numbers. In order to make the search for D € D efficient, it is better to 
control the class number beforehand. Tables can be made, but for larger computations, we need 
a fast way to compute h(—D). Subexponential methods exist, assuming the Generalized Riemann 
Hypothesis. From a practical point of view, our D’s are of medium size. Enumerating all forms 
costs O(h?) with a small constant, and Shanks’s baby-steps/giant-steps algorithm costs O( Wh) but 
with a large constant. It is better here to use the explicit formula of Louboutin [32] that yields a 
practical method in O(h) with a very small constant. 


5.2. An improved Cornacchia algorithm. Step 1 needs squareroots to be computed, then 
performs Cornacchia’s algorithm. Briefly, the algorithm runs as follows (see [38]): 


procedure CORNACCHIA(d, p, t) 
{t is such that t? = —d mod p, p/2 <t <p} 
8) 6. ph St ee = ee Se 
b) for i > 0 while r;_; > \/p do 
compute (a;;7;) Such that 9 Sage +7, 0S < Ties 
let Wi = Wi-2 + AW] («) ; 
c) if r?_, + dw?_, =p then return (r;_1, wj_-1) else return 9. 


We end the for loop once rj_1 < \/p < rj-2. As is well known, the a;’s are quite small and we 
can guess their size by monitoring the number of bits of the r;’s, thus limiting the number of long 
divisions. One can use a fast variant for this half gcd if needed, in a way reminiscent of Knuth. 

Moreover, we know that this algorithm almost always returns the empty set in step 2c), since 
the probability of success if 1/(2h(—d)). Therefore, when h is large, we can dispense of the multi- 
precision computations in equation (*). We replace it by single precision computations in base O° 
(b = 32 or b = 64): 


Wi = Wy_-2 + a;wj_1 mod og 


and at the end, we test whether r?_, + dw?_, = p mod 2°. If this is the case, then we recompute 
the w,’s and check again. 
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5.3. Factoring m. In practice, trying to factor m = 2N’, for a probable prime N’, is too restrictive; 
instead, we check whether m = cN’ for some B-smooth number c. The parameter B is critical. 
As shown in [20], the number of probable prime tests we will have to perform is heuristically 
t = O((log N)/(log B)) and we will end up with N’ such that N/N’ = log B. 
For small N, we can factor lots of m doing the following. In a first step, we compute 


ri = (N +1) mod p 
for all p; < B, which costs 7(B)L log B, where 7(B) = O(B/ log B) is the number of primes below 
B and the other term being the time needed to divide a multi-digit number by a single digit number. 
Then, sieving both m = N+1-—U and m’ = N+1+U is done by computing u; = U mod p; 
and comparing it to +r; for primes p; such that (—D/p;) #4 —1. See [3, 35] for more details and 


tricks. 
The cost of this algorithm for ¢ values of m is 


O(BL) + O(t- BL) 


where the second term comes from the computation of U mod p;, which is roughly twice as fast as 
that of (N + 1) mod p;, since U is O(v N). Since we need to perform also t probable prime tests 
(say, a plain Fermat one), then the cost is 


O(tBL) + O(tL?*“) = O( BL?) + O(L?**) 


and therefore the optimal value for B is B = O(L). 
When N is large, it is better to use the stripping factor algorithm in [20], for a cost of O(B(log B)?), 
the optimal value of B being B = O(L’). 


5.3.1. Remark. Suppose we have found N’ and that m= N+1-—U =cN’. Then we will have to 
compute 
r, = (N’ +1) mod p; = (r; — uj)/e + 1 mod 9. 


Computing the right hand side is faster, since ¢ is ordinarily small compared to N’. 


5.3.2. Using an early abort strategy. This idea is presented in [20]. We would like to go down as 
fast as possible. So why not impose N/N’ greater than some given bound? Candidates N’ need be 
tested for probable primality only if this bound is met. From what has been written above, we can 
insist on N/N’ & log B. Our implementation introduces a parameter 6 and imposes N/N‘ > 9°. 


5.3.3. Using new invariants. Tackling larger and larger N’s forces us to use larger and larger D’s, 
leading to larger and larger polynomials Hp. For this to be doable, new invariants had to be used, 
minimizing the height of their minimal polynomials. This task was done using Schertz’s formulation 
of Shimura’s reciprocity law [39], with the invariants of [18, 16] (alternatively see [25, 24]). Note 
that replacing 7 by other functions does not the change the complexity of the algorithm, though it 
is crucial in practice. 


5.3.4. Step 3 in practice. We already noted that this step is the theoretically dominating one in 
fastECPP, with a cost of O(L3+#+”). In practice, even for small values of h, we can assume v ~ 0 
(using for instance the algorithm of [33] for polynomial multiplication). 

One way to reduce the cost of this phase in practice is to use smooth values of h, and use Galois 
theory to factor Hp(X) over a tower of extensions. Then, we replace the factorization of a degree 
h polynomial by a list of smaller ones, the largest prime factor of h being logh. We already used 
that in ECPP, using [29, 17]. Typical values of h are now routinely in the 10000 zone. 

It could be argued that keeping only smooth class numbers is too restrictive. Note however, that 
class numbers tend to be smoother than ordinary numbers [10]. 
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5.4. fastECPP. We give here the expanded algorithm corresponding to step 1’. Using a smooth- 
ness bound B, we need approximately t = exp(—7) log N/log B values of m and therefore roughly 
t/2 discriminants. The probability that D is a splitting discriminant is g(—D)/h(—D). Therefore 


we build discriminants until 
S © g(—D)/h(—D) © t/2. 
D 


One way of building these discriminants is the following: we let r increase and build all or some 
of the subsets of {q},...,q} until the expected number of D’s is reached. After this, we sort the 
discriminants with respect to (h(—D)/g(—D), h(—D), D) and treat them in this order. 


6. BENCHMARKS 


First of all, it should be noted that ECPP is not a well defined algorithm, as long as one does 
not give the list of discriminants that are used, or the principles that generate them. 

Since the first phase of ECPP requires a tree search, testing on a single number does not reveal 
too much. Averaging on more than 20 numbers is a good idea. 

Our current implementation uses GMP [26] for the basic arithmetic, which enables one to use 
mpfr [28] and mpc [19], thus leading to a complete program that can compute polynomial H p’s on 
the fly, contrary to the author’s implementation of ECPP, prior to version 11.0.5. This turned out 
to be the key for the new-born program to compete with the old one. 

We give below some timings obtained with our new implementation, after a lot of trials. We 
used as prime candidates the first twenty primes with 1000, 1500, and 2000 decimal digits. Critical 
parameters are as follows: we used D < 10’, h < 1000, 5 = 12 (see section 5.3.2). For 1000 and 
1500 decimal digits (resp. 2000 digits), we limited the largest prime factor of h to be < 30 (resp. 
100). For the extraction of small prime factors (used in the algorithm described in [20] and denoted 
EXTRACT in the sequel), we used B = 8- 10°, 10", 3-107 for the three respective sizes. 

SQRT refers to the computation of the VE , CORN to Cornacchia, PRP to probable primality 
tests; HD is the time for computing polynomials Hp using the techniques described in [16], jmod 
the time to solve it modulo p; then 1st refers to the building phase (step 1), 2nd to the other ones; 
total is the total time, check the time to verify the certificate. Follow some data concering D, h and 
the size of the certificates (in kbytes). All timings are cumulated CPU time on an AMD Athlon 64 
3400+ running at 2.4GHz. In the tables, avg stands for average and std for standard deviation. 


Looking at the average total time, we see that it follows very closely the O((log N)*) prediction. 
According to our analysis, the main costs are the square roots, the pseudo primality test and finding 
a root of Hp, with the same O(L**“) heuristic cost. Square roots seem to become comparatively 
less important as D increases. 


7. CONCLUSIONS 


We have described in greater details the fast version of ECPP. We have demonstrated its ef- 
ficiency. As for ECPP, it is obvious that the computations can be distributed over a network of 
computers. We refer the reader to [20] for more details. Note that the current record of 15041 
decimal digits (with the number 44057638 + 2638" see [21]), was settled using this approach. 
Many more numbers were proven prime using either the monoprocessor version or the distributed 
one, most of them from the tables of numbers of the form «¥ + y* made by P. Leyland?. 


thttp://www.leyland.vispa.com/numth/primes/xyyx.htm 
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SQRT 
CORN 24 17 
EXTRACT 84 74 
PRP 124 102 
7 
99 
Ist 276 
2nd 136 
total 


nsteps 
certif 
D 


427 
140 
EXTRACT 282 
PRP 903 


171 

95 

230 

664 

3) 13 9 
868 1590 1192 
368 649 908 


1322 2240 1701 


nsteps 
certif 


‘TABLE 2. 1500 decimal digits 


Cheng [9] has suggested to use ECPP to help his improvement of the AKS algorithm, forcing 
m = cN’ to have N’ — 1 divisible by a given prime large prime of size O((log N)*). The same idea 
can be used to speed up the Jacobi sums algorithm, and this will be detailed elsewhere. 
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any 


aD oR wn 


10 


11 
12 


13 
14 
15 


16 


17 


SQRT 
CORN 
ee 


ae S57 4888 3778 


1398 2120 1777 


aa 4494 6795 5057 
Peak [a [aan] 2] a 
236 
1420 
9760387 | 285217 | 1026529 


‘TABLE 3. 2000 decimal digits 
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